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>- Abstract 
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A practical method for finding free energy barriers for transitions in high-dimensional classical and 
quantum systems is presented and used to calculate the dissociative sticking probability of Hi on a metal 
surface within transition state theory. The reversible work involved in shifting the system confined to a 
hyperplane from the reactant region towards products is evaluated directly. Quantum mechanical degrees 
Q\ ' of freedom are included by using Feynman Path Integrals with the hyperplane constraint applied to the 

centroid of the cyclic paths. An optimal dividing surface for the rate estimated by transition state theory 
is identified naturally in the course of the reversible work evaluation. The free energy barrier is determined 
relative to the reactant state directly so that an estimate of the transition rate can be obtained without 
requiring a solvable reference model for the transition state. The method has been applied to calculations 
of the sticking probability of a thermalized hydrogen gas on a Cu(llO) surface. The two hydrogen atoms 
and eight surface Cu atoms were included quantum mechanically and over two hundred atoms in the Cu 
crystal where included classically. The activation energy for adsorption and desorption was determined and 
found to be significantly lowered by tunneling at low temperature. The calculated values agree quite well 
with experimental estimates for adsorption and desorption. Dynamical corrections to the classical transition 



state theory rate estimate were evaluated and found to be small. 



1. Introduction 

The dissociative adsorption of molecules on surfaces of solids is of central importance in surface catalysis 
and has been extensively studied both experimentally and theoretically. Hydrogen adsorption on copper 
surfaces has become 'the classic example of activated dissociative adsorption of a molecule at a surface'. 1 
Many experimental studies have made use of molecular beams with molecules in a selected initial state 
impinging on the surface. Theoretical studies have focused on classical trajectory calculations and quantum 
wavepacket propagation to explore the dissociative sticking dynamics given a well defined initial state of 
the molecule. 1 A great deal of qualitative insight has been gained, but many questions remain open. On 
the theoretical side, the problem is that full quantum mechanical treatment of all six hydrogen degrees of 
freedom is impractical with present computers and currently available techniques for quantum wavepacket 
propagation. The theoretical work has therefore largely been confined to quantum calculations on lower 
dimensional model systems 2-4 or mixed quantum and classical treatment. 5 Not only is a six dimensional 
wavepacket propagation a big order, but such a calculation would need to be repeated many times to 
average over the surface vibrational degrees of freedom to get a single value of the sticking probability. Low 
dimensional simulations have demonstrated clear quantum effects, 1,2 and calculations in higher dimensions 
have indicated that dimensional effects might be equally important. 3,4 

We report here on a fully quantum and thermally averaged simulation of the dissociative sticking of 
Hi on Cw(llO). The large number of degrees of freedom can be included at the expense of the amount of 
information obtained about the dynamics. We apply transition state theory, a statistical theory of rates, 
to this problem and study entropic and quantum effects on the sticking probability. Quantum effects are 
included by Feynman path integrals (FPIs) which implicitly include thermal averaging over quantum states. 6 
The simulation mimics a thermodynamic experiment, where the impinging molecules have a thermal energy 
distribution at the temperature of the substrate. This experiment has been carried out in the laboratory 
of Campbell and coworkers 7 who obtained the thermodynamic activation barrier for dissociative adsorption 
from the temperature dependence of the sticking probability. Our simulations using a rather simple model 
potential surface show a clear onset of a quantum mechanical regime as the activation energy drops at a 
temperature around 400^. The experimental measurements of the sticking coefficient were taken slightly 
above this transition temperature. 

The method we present here involves evaluation of reversible work to determine the free energy of the 
system along the reaction coordinate and, in particular, the free energy barrier for the transition. Traditional 
methods for evaluating free energy differences between two states, in particular a transition state and a 
reactant state involve integration (using Monte Carlo sampling) over a parameter A which smoothly converts 
the potential from one state to another. 8 The advantage of the present technique is that the transition state 
need not be known beforehand, and since the integration is over the reaction coordinate the intermediate 
states in the calculation have some physical significance. Voter 9 has presented a technique for direct Monte 
Carlo calculation of free energy differences between separated, localized states by translating one to overlap 
the other. This method eliminates the need for sampling intermediate states but again requires a priori 
knowledge of the transition state, and can have sampling problems if the transition state potential is shaped 
much differently than a slice of the reactant state potential. 

The method described here has previously been applied to calculate transition rates in a model system 
consisting of an Eckart barrier coupled to a harmonic oscillator 10 . Even in such a small system, the reversible 
work evaluation was more efficient than calculation of the partition function with a harmonic reference. 
Here we give expressions directly applicable to large systems and apply the method to hydrogen dissociative 
adsorption. Preliminary results of these calculations have been presented elsewhere. 11 

The organization of this paper is as follows. We first discuss classical Transition State Theory (TST) in 
section 2. Then quantum TST is discussed in section 3. In section 4 we present a derivation of the equations 
required for calculating free energy barriers in classical systems by evaluation of reversible work. The 
application of these equations involves statistical sampling of forces acting on the system while constraints 
are applied. Following Gillan, 12 we then generalize these expressions to quantum systems in section 5 by 
representing quantum particles by Feynman path integrals and applying the constraints to the path integral 
centroids. In section 6 we discuss the hydrogen dissociative sticking problem and present our results on the 
Hi — Cu(110). A summary is presented in section 7. 
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2. Classical Transition State Theory 

Transition state theory is a method for estimating rate constants of transitions, for example, chemical 
reactions and diffusion events. 13 ~ 15 It is applicable to slow transitions or rare events. It is a statistical theory, 
but is based on dynamics. The underlying idea is that the equilibrium density at the reaction's transition 
state - a statistical quantity - is proportional to the reaction rate - a dynamical quantity. A central aspect 
of transition state theory is the identification of a dividing surface separating the initial state (the reactants) 
and the final state (the products). Such a dividing surface is frequently referred to as the 'transition state'. 
All possible transitions from the initial state to the final state take the system through the dividing surface. 

A well-chosen dividing surface also separates the initial and final states in terms of dynamics. A system 
at such a transition state may become either reactants or products, and once it becomes one or the other, 
it is likely to remain so for a much longer time than it spent near the transition state. Thus, the system's 
character in the future is strongly dependent on what it does during its brief visit to the transition state. 
This makes the transition state and what happens there of central importance. 

The system consists of the degrees of freedom of primary interest and all degrees of freedom to which 
they are coupled. For instance, in the H2 dissociation on Cu, the system consists of both the dimer and the 
Cu lattice. The hydrogens move furthest and are the primary degrees of freedom, and the reaction's progress 
- the reaction coordinate - is determined by the hydrogens' positions. The surface degrees of freedom are 
also specifically included but do not enter the reaction coordinate or the definition of the transition state; 
they are effectively a bath. 



2.1 The Central Assumptions of TST 

The first central assumption of transition state theory is that the reactants have a Boltzmann distribution 
of energy in every degree of freedom and that the rate of transition is low enough that this distribution is 
maintained. 

As a consequence, the system will have a Boltzmann distribution of energy at the dividing surface, 
even if the system doesn't exchange energy with a bath during its climb up the barrier. The density at the 
barrier is low because the barrier turns most trajectories around. Those trajectories which do make it to 
the dividing surface, however, have a momentum and spatial distribution at the barrier which is Boltzmann 
at the reactant temperature except that states originating from products are missing. 1 ' It is important 
to reallize that TST does not need to make any assumptions about the ability of the system to achieve an 
equilibrium distribution at the dividing surface. 

If the reactants are specifically prepared in a non-Boltzmann distribution, for example, in a molecular- 
beam experiment where different reactant modes may be at different temperature, or the surface and adsor- 
bate temperatures are different, then the system - consisting of both the incident molecule and the surface 
- is not at thermal equilibrium, and TST does not apply. We calculate here the sticking probability of a 
thermalized gas, where the gas and surface temperature are the same. 

The second central assumption of transition state theory is that the system crosses the dividing surface 
only once in going between the reactant and product states. A reactive trajectory is one which takes 
the system from steady-state reactants to steady-state products. The actual reaction rate depends on the 
density of reactive trajectories at the dividing surface. Transition state theory assumes that any crossing of 
the dividing surface in the direction from reactants to products is a reactive trajectory. In reality, trajectories 
can cross the dividing surface many times. Any trajectoriy that crosses the dividing surface an even number 
of times, starting from either side, is not reactive, but is counted in TST. Further, a trajectory which crosses 
several times before ultimately reacting is reactive, but is counted several times in TST while only forming 
one product. As a result, the TST estimate of the rate is always an overestimate, but in favorable cases and 
for a good choice of the dividing surface TST can give an estimate that is nearly exact. 

The fraction of trajectories that have been overcounted in the TST estimate can be found only by 
studying the dynamics of the system. Transition state theory assumes that the probability, k, of continuing 
to products, once the dividing surface is reached with a momentum directed towards products, is unity. If 
the system sometimes recrosses the dividing surface, k < 1. For a given choice of the dividing surface a 
transmission coefficient, k, can be found such that the exact rate constant is k exact — n k T . 



For a good choice of the dividing surface, k is as large as possible, and k TST is as close as possible to 
the actual rate constant. A dividing surface far into the reactant state would result in many crossings, most 
of which would lack enough energy to cross the barrier and would soon recross it. For this dividing surface, 
k<1, and k T is a poor estimate of the rate. Similarly, a dividing surface far into the product state would 
be a poor choice, because most of the density there comes not from paths that have recently reacted, but 
from products that have recently tried and failed to back-react. 

While TST predicts that the sticking coefficient depends only on interatomic forces and not on atomic 
masses, the recrossings can be mass dependent, for example, if the products re-desorb because energy is 
not redistributed between different modes on the time-scale of the reaction. This phonon coupling is mass- 
dependent, and hence can lead to kinetic isotope effects that are not predicted by TST. 

Below, we describe studies we have carried out of dynamical trajectories for Hi dissociative sticking 
on Cit(llO). We found that dynamical corrections play a minor role in this system. TST can therefore be 
expected to give reliable estimates of the sticking coefficient. 



2.2 Variational Transition State Theory 

Since classical TST gives an upper bound to the rate, and can be applied to any dividing surface between 
reactants and products, it is possible to search for the dividing surface for which the estimated rate becomes 
a minimum; this estimate is closest to the true rate. This dividing surface will have the fewest recrossings. 

This optimization of the TST estimate is called Variational Transition State Theory (VTST) 17 . The 
variation of the plane can be done separately at different temperatures, or at different initial reactant 
momenta, or by making the dividing surface non-planar, any added flexibility giving an improved rate 
estimate. There exist hypersurfaces that separate the dynamics completely; in general, these are curved, 
and different for each reactant energy. 13 ' 18 If such a hypersurface could be found, the rate associated with 
it would be the exact reaction rate: there would be no recrossings. However, for high-dimensional systems, 
finding an optimal surface, even within a restricted set (such as hypcrplanar surfaces), is a non-trivial 
problem. 

2.3 The Transition State Partition Function 

Figure 1 shows a reaction path, T(s), between reactants and products, parameterized with reaction 
coordinate s. If N p is the number of atoms of primary interest; that is, those that undergo a significant 
displacement during the transition, then T s is a 37Vp-dimcnsional vector pointing at position s on the path. 
The reaction path is defined in terms of these coordinates and is independent of all others. The dividing 
surface, Z+, is a (3N p — 1) -dimensional cut through the 3iVp-dimensional system. In the case of the H2/CU 
system, only the hydrogens change position significantly during the reaction, so the dividing surface is a 
5-dimensional cut through the 6-dimensional Hi system. All degrees of freedom independent of T comprise 
the bath. In the -fl^/Cu system, the bath consists of all the Cu atoms. If, for instance, two of the Cu atoms 
moved significantly to enable the dissociation to happen, the reaction path would include those Cu atoms, 
and all other Cu atoms would be the bath. 

The dividing surface has one dimension fewer than does the reactant state. A transition state with the 
same dimensionality as the reactant can be defined by arbitrarily introducing a small width S to the dividing 
surface. Eventually this width is cancelled out and does not enter the final expressions. The probability of 
finding the system at the transition state is Q* /Q R . The transition state partition function can be written 
as Q* = Q^5, where Q^ is the partition function, taken over the dividing surface, over all system coordinates 
other than the reaction coordinate. 

If the dividing surface, Z?, is specified to be a hyperplane (see fig. 1) with unit normal vector n 
intersecting the path at r+ and letting r be the system coordinate in the hyperplane, then the partition 
function is 

<oi= f e -mr) 5 



n-(r-r+) dv (1) 

where V is the interaction potential, /? = l/fceT, Ub is the Boltzmann constant and T is the temperature of 
the system. The Dirac 5— function imposes the hyperplane constraint. We will deal only with hyperplanar 
dividing surfaces in this paper. 



2-4 The Reaction Rate 

Once a system has reached the transition region, its rate of escape is 1/t = v±/<5, where t is the time 
it spends there; vj_ = v • n is its velocity perpendicular to the hyperplane. To find the forward flux, only 
those points in the transition region that are forward-bound should be counted. These comprise half the 
population there. The escape rate constant is Q* /Q R times the average of 1/t: 

The factors of 5 in 1/t and Q+ have cancelled. Only positive velocities count because back reactions are 
ignored. The average forward velocity is 

(Kl) = k B T 
2 V2tt ^ksT 

The TST estimate for the rate constant is 



>TST __ k BT Q+ 



y/2TT^k B T Q R ' 

Here /i is an effective mass for the reaction coordinate. 



(4) 



3. Quantum Transition State Theory 

Several different extensions of classical TST to quantum systems have been proposed. 12 ' 19-27 We will 
focus here on methods where the classical statistical averaging is replaced by quantum statistical averaging. 
3.1 Quantum Statistical Mechanics 

Fcynman and Hibbs 8 describe a reformulation of quantum mechanics in terms of Feynman path integrals 
(FPI). This method involves integrating over all possible paths a system could take from one place to another. 
If done in real time, it gives quantum dynamics. In imaginary time, it gives quantum statistics. In complex 
time, it gives a combination of these. 

The quantum mechanical partition function, Q, can be written as 

Q= [ e - v °ffM T W kBT Dx(T), (5) 



where Dx{t) includes P integrals over coordinates x(t) ~ {xi,X2, ■ ■ ■ ,Xp} comprising cyclic paths. The 
paths are represented to resolution P in this discrete approximation. Each path contributes to Q according 
to its effective potential energy, V e ff 8:28 



v.tt - E 
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i=l 



P 



(6) 



with the spring constant given by k spr = mP/h [3 where m is the particle mass and f3 = l/ksT . Eqn. (5) 
becomes exact as P — » oo. 

The classical system has effectively been replaced by P images of itself. Each image is subject to forces 
1/P as strong as the classical system and interacts with its neighbors through the harmonic springs. In 
the classical limit as the mass or temperature become larger, k spr — > oo, and the stiff springs prevent any 
derealization. The P images combine to give back the original classical system. 

With P = 1 the classical approximation is recovered. For P > 1, the springs allow the chain to 
delocalize, which gives quantum statistical effects. If the chain is at a barrier, it can drape down on cither 
side to sample a lower effective potential; this represents tunnelling. If the chain is confined in a small area 
(as is a particle in a box), it will have fewer available arrangements, which gives it a lower entropy. This 
confinement entropy corresponds to zero-point energy. 



3.2 Centroid Density based QTST 

A rate theory for transitions in quantum systems based on FPI was introduced by Gillan. 12 He calculated 
the thermal equilibrium probability density of the quantum particle at various locations by introducing a 
constraint on the system, fixing the centroid of the FPI chain 

I p 

xo = — Y^ Xi ( 7 ) 

i 

at a given point, and taking the average over the remaining quantum degrees of freedom and classical bath 
degrees of freedom. In particular, he evaluated the probability of finding the centroid at the saddle point of 
the potential surface between reactants and products, relative to the probability of finding the centroid in 
the reactant region. By gradually shifting the centroid constraint from the reactant region up to the saddle 
point and monitoring the reversible work, he obtained a free energy difference. A rate estimate was obtained 
by assuming the rate is proportional to the relative density of the centroid at the saddle point. 

In one dimensional systems the classical limit of this estimate agrees with classical TST. However, in 
higher dimensions it does not. The centroid in Gillan's formulation is confined to a point rather than to 
a dividing surface; the transition rate is related to the centroid density at two specified points. 10 ' 12 This 
estimate includes bath degrees of freedom, and various FPI chain configurations of the quantum particle 
(thus it picks up some tunnelling and zero-point energy effects), but it does not include averaging over 
centroid degrees of freedom perpendicular to the reaction path. Thus, it does not fully include effects due to 
varying width of the reaction channel. 

Following Gillan's idea, Voth, Chandler, and Miller (VCM) 22 provided a quantum generalization of 
traditional TST by replacing the classical coordinate in eqn. (1) by the centroid coordinate of the quantum 
particle. Again, a constrained FPI formalism was used, but with the centroid restricted to a multidimensional 
dividing surface rather than to a point. The calculations were done by explicit evaluation of the partition 
function, using an analytically solvable reference problem. This method was applied to calculations on a 
restricted geometry of the reaction H + H2 — * H2 + H and found to give results in good agreement with 
more accurate calculations. 

Messina, Schenter, and Garrett (MSG) 24 ' 25 generalized the VCM expression and allowed for variational 
optimization of non-planar and momentum-dependent dividing surfaces. The dividing surface was chosen to 
minimize the calculated rate obtained by direct evaluation of the dividing surface partition function. MSG 
is a quantum analogue of classical VTST. In these calculations a Monte Carlo sampling procedure is used 
to find the partition functions explicitly. 

These centroid density methods have been applied to various test problems and have been shown to give 
good rate estimates. Transitions in three dimensional systems involving one quantum particle have also been 
studied using these techniques. Gillan studied the diffusion of a hydrogen atom in bulk metals. 12 Mattsson, 
Engberg and Wahnstrom 27 and Sun and Voth 28 studied the diffusion of a hydrogen atom on Ni(100) and 
Cu(100) surfaces, respectively. 

We describe below a different computational method, reversible work based TST, for carrying out 
classical and quantum TST calculations. For classical systems the results are equivalent to VTST. For 
quantum systems, the results are equivalent to MSG. The advantage of the reversible work formulation is 
that the method can more easily be applied to high dimensional systems. After describing the method and 
its foundation, first as it applies to classical systems (section 4) and then as it applies to quantum systems 
(section 5), we then describe an application to Hi dissociative adsorption, a reaction path involving two 
quantum particles (the H atoms) and a bath involving eight quantum particles (surface Cu atoms) and a 
couple of hundred classical particles (section 6). This is by far the largest application of QTST we are aware 
of. Our conclusions arc given in section 7. 



4. Reversible Work Formulation of classical TST 

k TST is written in terms of a partition function ratio, or a free energy difference. The free energy can 
be calculated by the method of reversible work rather than by calculating the partition functions explicitly. 
The rate constant, eqn (4), can also be written 



,tst = k B T Q zR Q% 

y/2nfik B T Q R Q ZR ( ' 

where we have introduced Q z a partition function for the system confined to a reference hyperplane Z R 

in the reactant region. Q z /Q R has units of inverse length. ~~ zR is unitless and, if there are no recrossings 

of the transition state, is the reaction probability: (number of reactions) /(number of attempts). The rate 
is (average positive velocity) x (average reciprocal length) x (reaction probability). For surface adsorption, 

-^tt is the sticking coefficient. 



4-1 Relative probability and free energy 

The reaction probability can also be written in terms of the free energy,yl, of the transition and reference 
planes: 

Q$/Q zR =e-^-* zR \ (Q) 

AA = A* - A z " 

In the reversible work evaluation the hyperplane is moved stepwise from Z B to Z+, as indicated in figure 2, 
and the free energy difference, AA, is accumulated at each step. 

Figure lb shows a progression of hyperplanes Z s defined along the entire reaction path T s between 
reactants and products. Given a path, the sequence of hyperplanes can be constructed by choosing the 
tangent vector of the path ^T s , at various points along the path as the normal to a hyperplane associated 
with that reaction coordinate 



-Ts/ 
ds 



ds 



(10) 



The progression of the hyperplane has to be gradual enough that a first order evaluation of the free energy 
difference will be accurate. As the reaction coordinate is incremented by ds, where 

ds =\\T s+ds -r s \\, (ii) 

the hyperplane intersection with the path moves from T s to r s +d s and the normal vector changes from n s 
to n s+ ds- 

It is convenient to introduce a coordinate system aligned with the hyperplane. These local coordinates, 
z, are obtained from the actual coordinates, r, by applying a unitary rotation matrix U: 

r s - r s = u s z, (i2) 

Ys+ds — r s+ rf s = U s +dsZ- 

We choose the first unit vector in the local coordinate system to be normal to the hyperplane. The hyperplane 
constraint, (r s — F s ) ■ n s = 0, can then be written as z\ =0. The first column of U s is n s . 

The free energy difference corresponding to an increment ds in the reaction coordinate is dA — — k B T ln(Q s+( j s /Q s )| 
where 

r e- 0v -^6{zi)dz 

(13) 
Qs+ds = I e-* v '+"W5(*i)dz. 



From the point of view of a system confined within the hyperplane, the potential energy changes as the 
reaction coordinate is incremented 

V s (z) = V(r s ) 

V s+ds (z) = V(r s+ds ). 
6 



This ramping of the potential energy of the system is illustrated in fig. 2. Since the planes are nearby, 

V s+ds (z) = V s (z)+5V(z), (15) 

where SV(z) is small for all values of z the system is likely to visit. 
The partition function ratio is 



which can be written as 



W s-\-ds 



Qs 



Q s fe-P v ^)S( Zl )dz 

-f 3Sv ^) s =e-^ sv ^°+0(5V 2 ) 



(e 



(16) 



(17) 



since SV is small. The average is calculated in plane Z s . Thus, if the partition functions are written in terms 
of the constant internal coordinate z, moving the hyperplane from one position to another is, in effect, the 
same as ramping the potential energy seen by a system confined to the hyperplane by an amount 5V(z). 
The free energy change becomes simply the thermal average of the local internal potential ramp SV(z): 



dA=(5V(z)) s +0(SV 2 ). 
The external coordinates of the system in the two hyperplane locations are related by 

T s +ds = U s+ d s U s (r s — T B ). 



r s +ds 



(18) 



(19) 



Writing U s + ds = U s + 5U S , U s + ds Uj = U s Uj + 8U s Uj = I + SUgUj' and substituting into the equation 
above gives, after rearranging, 



r s+ds - r s = (T s+ds - T s ) + 5U s Uj(r s - T.) 
= n s ds + 6U s Uj(r s - r s ) 



(20) 



using the fact that we chose the normal to be tangent to the path, eqn (10). The first column of U s is n s . 
The other columns can be assigned arbitrarily; the free energy of the plane is independent of them. These 
column vectors Uj comprise an orthonormal basis for the hyperplane. The geometrical difference between 
the two locations of the hyperplane is entirely specified by T s , T s+ds , n s , and n s+d s- We will now show that 
the free energy change can be expressed entirely in terms of this geometry, without explicit calculation of 
the internal basis vectors. 

If D is the dimensionality of the primary system (3N p in a three dimensional system), then 



6U S = [dn s du.2 



dur 



(21) 



Since u^ • Uj = 5ij for an orthonormal basis, ^j(u^ ■ Uj) = 0. This property makes the matrix 5U a U] 



antisymmetric. 



SUsVZ 



{ dn s du 2 



dup] x 



u 2 



(22) 



Writing, for convenience, n s as ui, and making use of the antisymmetry, this can be written as 



6U s Uj = 



- J2 du i2 ua . . . 

-J2 du ^ u t2 

-Y^dUnUiD -Y,dUi2UiD ■■■ 
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- J2 du iD u iX 

- J2 du iD u i2 



(23) 



where Uij is the jth component of the zth basis vector. 

Operating with this matrix on the vector R = r s — T s yields, after rearrangement, 

6U.U?(R) = -Y u t (du t • R). (24) 



i 

The internal potential energy ramp is 

SV(z) = V(r s+ds ) - V(r„) = -F(r s ) • (r s+ds - r s ) 

= -F(r s ) • n s ds - F(r s ) • (SU s Uj(R)) 

(25) 
= -F(r s ) • n s (ds - dn s • R) + J^( F ' u *)(^ • R). 

The thermal average of each term in the sum over the in-plane basis vectors is zero, ((F-Uj)(dUj-R)) = 0. 
Since the unit basis vector u^ is orthogonal to its change dvn 1 the average of the product is the product of 
the averages. ((F • Uj)) = for any in-plane u^, whether or not any orthogonal coordinates are specified. 
(dui • R)(F • Uj) is therefore individually zero for each specified (dui ■ R), and the overall average vanishes. 

The angle between the planes' normal vectors, n s and n s+c is, is d6 = ||rfn s ||. The first-order curvature 
of the reaction path r at s is k — d6/ds. We write for the normal component of the force, F n = F(r s ) • n s , 
and for the turning-direction component of the distance from the reference point, Rt = (r s — T s ) ■ dn s /dd. 
Then we can write the free energy difference between the planes as 

dA(s) = - (F n (l-nRt)), ds . (26) 

Here —(F n R t ) s is the average torque opposing the rotation of the plane around the reference point T s in the 
direction dn s /d9. The amount of work required to rotate the plane by d9 in this direction is 

dA rot (s) = - (F n R t ) s d0 = - K(F n R t ) s ds. (27) 

The free energy difference between the reference hyperplane, Z R , and a hypcrplanc Z s corresponding 
to reaction coordinate s is 

AA(s) = - ( (F n (l- nR t )) s , ds 1 . (28) 

JR 

If all coordinates are specified, the instantaneous force is the mean force. The line integral of the 
instantaneous force is the potential energy. The reversible work involved in shifting the system from one 
completely specified configuration to another is just the potential energy difference. 

4-2 Implementation of the method 

First a reaction path leading from the reactant to the product region is defined. A sequence of hyper- 
planes intersecting the path is then defined. The hypcrplanes must be spaced closely enough for the first 
order expression, eqn. (26), to be accurate enough. The average force and torque acting on the system 
confined to the hyperplane are then statistically sampled and the results integrated according to eqn. (28). 

A convenient choice for the path is the minimum energy path (MEP), but any path that does not lead to 
large numerical cancellations in the force evaluation can be used. The hyperplane normals can conveniently 
be chosen to be tangent to the path (as we have done here), but this is not necessary. To collect the average 
force acting on the system in each of the hypcrplanes, a Monte Carlo or a molecular dynamics simulation 
can be carried out subject to the constraint that the system remains in the hyperplane. More specifically, 
considering a primary system consisting of N p atoms in 3-dimcnsions, a SA'p-dimensional unit vector, the 
hyperplane normal n s , is used to constrain the dynamics. At each time step the component of the force 
vector normal to the hyperplane, F • n s , is collected for the evaluation of the thermal average, < F n >. But 
then, this component is projected out of the force to give a new, revised force F rev — F — F • n s which is 
used in the simulation giving the statistical sampling. 

We have used the velocity Verlet algorithm to carry out the molecular dynamics sampling of canonical 
phase space. 31 The constrained dynamics conserve energy as long as the constraint vector is constant. The 



optimal dividing surface Z^ for the TST rate estimate is the hyperplane Z s corresponding to the maximum 
free energy, where force < F n > s and the torque < F n R t > s acting on the hyperplane both vanish. It could 
be necessary to adjust the orientation of the hyperplane at and near the transition state in order to zero the 
torque. Alternatively, the path can be redefined to construct a 'zero torque path' as we have done for an 
Eckart barrier test problem. 10 As wc will show below, the torque acting on the MEP-constructed hyperplane 
in the hydrogen dissociative adsorption simulation turns out to be very small, leading to negligible rotational 
contribution to the free energy barrier. 

Each plane has a certain free energy, AA(s), relative to the reference plane in the reactant state. This 
free energy gives a TST estimate for the rate, treating that plane as the dividing surface: 

i,tst _ ^bT Q s _ k B T Q pAAjs) / 9 q\ 
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The plane for which this estimated rate is lowest is the planar dividing surface at which there are the 
fewest recrossings. So VTST is accomplished automatically by the reversible work calculation. At some 
plane crossing the reaction path, the free energy will be a maximum; this gives the VTST rate estimate. 
The automatic inclusion of VTST is an advantage of the reversible work method. 

The free energy barrier found this way is path-independent as long as the planes near the barrier are 
good dividing surfaces. The optimal planar dividing surface is the one corresponding to maximum free 
energy, and therefore has zero average force and torque acting on it. If there is a net torque on the best 
plane of the sequence, it can be rotated against the torque, picking up some additional free energy, until the 
torque is zero. 

The idea of using a one-parameter sequence of generalized dividing surfaces and minimizing the cal- 
culated rate constant with respect to the parameter has been implemented by Garrett and Truhlar in a 
semiclassical version of variational transition state theory. 17 Although their methods employ different calcu- 
lations of the rate than the current work, the idea of optimizing the rate along a path is equivalent to the 
determination of the hyperplane giving the maximum free energy in the present method. 



5. Reversible work formulation of quantum TST 

The formulas derived for classical TST will now be generalized to quantum systems by using quantum 
statistical mechanics instead of classical statistical mechanics to estimate the free energy barrier. The central 
question is how to incorporate the constraints on the system; i.e. how to confine the system to a hyperplane. 
Following Gillan 12 wc represent each quantum particle in the system with a Feynman Path Integral and 
apply the constraints to the centroids of the path integrals instead of the atom coordinates as in the classical 
case. 

The hyperplane constraint then becomes 

Z s = (f - r a ) • n s = 0, (30) 

where ?o is a vector containing the controid coordinates of quantum atom FPIs. Using a discrete represen- 
tation of the paths, the centroid for atom k is 



i p 



- fc « (31) 



where the sum extends over all the FPI images of atom fc. The reaction's progress is written in terms of the 
FPI centroids of the quantum mechanical atoms in the primary system. Therefore, the same reaction path 
can be used as in the classical calculation. 

The quantum mechanical partition function for a hyperplane at s is: 

Qs = [e- 0v °"6[n s -(r o -r s )]Dr(T). (32) 



with V e ff given by eqn. (6). Just as in the classical case, the change in the hypcrplane free energy AA(s) 
can be evaluated. The quantum mechanical generalization of cq. (28) is 



Ai(s) = ~ f { otl F "( J ) (1 - «#*(*))>«' da' (33) 
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where the hyperplane statistical average is 
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(34) 



In the Appendix we prove that the in-plane components of the force on the images do not contribute to 
AA(s), just as in the classical case. 

The integration along the reaction path gives a free energy maximum and thereby the free energy barrier 
for the transition in the quantum system. The transition rate is then obtained from the free energy barrier as 
in the classical case. Again, the reversible work calculation automatically includes variational maximization 
of the hypersurface free energy. The transition rate obtained by this reversible work method is exactly the 
same as the variationally optimized rate for a hyperplanar dividing surface obtained in the MSG theory, 2 ' 25 
but the advantage of the present method is that a solvable reference system for the transition state is not 
required. 

There is no rigorous justification for this formulation of QTST. Other formulations of a quantum gen- 
eralization of TST have been proposed. The relationship between the estimate obtained from the QTST 
presented here and the exact reactive flux 32 is unclear. The rate constant for a quantum system calculated 
by this and other centroid methods has not been proven to give an upper bound to the rate in general. It 
has been found that optimization of the dividing surface leads to a strict upper bound on the rate in the 
harmonic limit as well as in a symmetric Eckart barrier test problem. 10,24 However, kQ T has in some other 
cases been found smaller (but only slightly smaller) than k exac t- 



6. Dissociative Sticking of Hydrogen on Surfaces 

We now turn to the application of the reversible work TST method to the dissociative sticking problem. 
The calculation of the sticking coefficient described here include all six H2 degrees of freedom quantum 
mechanically as well as eight Cu surface atoms and full thermal averaging over several hundred substrate 
degrees of freedom. However, the calculations are done within TST and therefore apply to sticking of 
thermalizcd hydrogen gas, not molecular beam experiments with selected initial states of the molecules. The 
large number of degrees of freedom can be included at the expense of detailed dynamical information about 
the system. 

Campbell, Domagala, and Campbell 9 measured the sticking coefficient of thermally equilibrated hydro- 
gen gas on Cu(llO) experimentally. These measurements give an activation energy of 0.62 eV in the range 
500 K to 700 K. Molecular beam experiments convolved to yield a thermal distribution give similar results. 33 
Also, the activation energy for thermal desorption of H2 from Cu(llO) has been measured by Anger, Winkler, 
and Rcndulic. 34 The results of our simulation are compared to these experimental measurements below. 



6. 1 The 'hole model ' for sticking 

While TST is most often applied to transitions involving large energy barriers, the foundation of TST 
requires only that the rate of transition is slow enough to allow the remaining reactants to stay at thermal 
equilibrium; and that once the dividing surface has been crossed, it will not be recrossed on the time scale 
of the transition. A bottleneck in phase space involving only an entropic barrier can satisfy these conditions 
- an energy barrier is not necessary. 

As an example, the TST assumptions apply well to an ideal gas in a box effusing through a small hole. 
There is no energy barrier; just an entropy barrier. Most atoms colliding with the wall will fail to escape 
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because the hole is much smaller than the box's cross-section. If the box volume is V and the hole has cross 
sectional area A, then Q+ /Q R = A/V and the TST estimate for the effusion rate is 

Nk TST = Apk B T = AP 

\J2-K\xkBT \J2-K\xkBT 

which is precisely the rate given by kinetic theory. 

A successful model of molecular adsorption, the 'hole model', 35 ' 36 is analogous to the effusion problem. 
The assumption there is that a molecule approaching the surface with the right orientation will adsorb. 
However, most orientations of the molecule are unfavorable and lead to a large potential barrier which 
reflects the molecule back into the gas phase. The transition rate is then estimated as the ratio of the 
volume of the hole relative to the total volume available to the incoming molecule. If the hole model is a 
reasonable model for dissociative sticking, then TST can be expected to give a good estimate of the sticking 
probability. 
6.2 The H 2 — Cw(llO) Potential Energy Surface 

A simple form for the potential function was used to speed up the calculation. An EAM-type 37 potential 
form was chosen and fitted to data for the Hi dimer and the chemisorbed hydrogen. It was also fitted to 
a grid of points taken from a semi-empirical LEPS potential 38 , which is expected to have the right shape 
qualitatively, as it includes overlap and exchange effects. However, the potential barrier was adjusted to 
have a much higher saddle point than in ref. 38. The potential barrier was scaled up to 0.688 eV, close to 
the estimate of Hand and Harris based on ab initio calculations. 39 (See table 1). 

The EAM form works well for copper. With a filled d band, Cu exhibits little directional bonding. Most 
of the electron density seen by a neighbor is due to the 4s electron. 

The EAM potential form is: 

f = X>(r«) + !>(/*). (36) 

ij i 

where is a pair potential that represents the screened coulomb repulsion between the ion cores, including 
the inner electrons, and F(p) is the energy of embedding an ion core into a valence electron gas of density p. 
The electron density pi at atom i is the sum of the contributions of all its neighbors, j: 

p* = Y.p a ^)' ( 3? ) 

3 

where p A (rij) is fitted to be the Hartree-Fock atomic electron density. 
We have chosen the following functional forms 

0(r«) - D A e- aAr ^ + D B e- aBr *i 

M 

F(Pi) = E f™P? ( 38 ) 

m— 1 

p A ( rij ) = Sr n (e-P^n + 7 e- /3sr «) 

4>(rij) is a pair potential that represents the screened coulomb repulsion between the ion cores. It is 
parametrized as a double exponential. Each type of pair interaction (Cu-Cu, H-H, and H-Cu) has its 
own parameter set {Da, a a, Db, cub}- The parameters for the pair potential, the embedding energy and the 
electron density are given in table 2. 

The pair potential and the electron density are shifted in such a way that they become zero at and 
beyond a cutoff distance, R cu t, chosen to be 6.lA. 



6.3 Finding Minimum Energy Path using Nudged Elastic Band method: 

The MEP is defined by the boundary conditions, and by the requirement that the force at any point be 
parallel to it. That is, for any direction p s perpendicular to the path, p s • -§^T S = 0, the potential energy is 
at a minimum: p s • VV^(r s ) = 0. 

If 



The boundary conditions are that the MEP must start at the lowest-energy reactant state and finish in 
the lowest-energy product state. For the H2/CU problem, the former has the dimer distant from the surface 
at its equilibrium separation; the latter has the hydrogens on neighboring sites in the same groove on the 
(110) surface. 

The 'Nudged Elastic Band' method 40 was used to find the path. This method finds the MEP varia- 
tionally. An elastic band, made of beads connected by harmonic springs, is strung between the reactant 
and product states. Using a set of images or replicas of the system to define a discrete transition path, the 
problem is turned into a minimization problem by defining an object function 

p- 1 p kP 

F(n,r 2 ,...,rp_i) = J2 V ^ + E — ( r *-r*-i) 2 (39) 

where the sum is over the 'true' potential of all the intermediate images of the system, and the second sum 
is an 'elastic glue' or 'spring energy' that keeps adjacent images together The end point images r and vp 
are fixed. Although this chain is in some ways analogous to the off-diagonal density matrix FPI, the MEP is 
found with no reference to quantum mechanics. Here, the spring constant is arbitrary except that it should 
scale as kP to ensure convergence to the MEP as more images are introduced in the chain, P — > oo. 27 

In order to accelarate the convergence of the chain to the MEP we introduce 'nudging'. The component 
of the net spring force on a bead parallel to the local path tend to keep the bead equidistant from its 
neighbors; whereas the perpendicular spring force tends to make it collinear with its neighbors. When the 
MEP is curved, the chain will tend to pull the chain off the MEP by the perpendicular spring force (thereby 
inducing corner cutting). We therefore zero the perpendicular component of the spring force. Each bead 
also feels a force due to the gradient of the 'real' potential energy. The component of this gradient force 
perpendicular to the path tries to pull the beads down toward the MEP; whereas its parallel component tries 
to pull the beads down toward either the reactant or product state, away from the saddle point. This effect 
would result in reduced resolution at the saddle point. We therefore zero the parallel component of this 
gradient force. Moving each bead in the direction of its adjusted force until the adjusted force becomes zero 
will put the beads on the MEP (to within the chain's resolution), equally spaced in the 3A p -dimensional 
space. The technique is not a global search: if there is more than one MEP, only one will be found. 

Figure 4 shows a MEP for the H2/CU problem. The initial state is asymmetric with the final state. Even 
so, the path becomes symmetric in the vicinity of the surface, lining up with the groove on the (110) surface. 
The Nudged Elastic Band method has therefore shown that the MEP is symmetric. In a six dimensional 
system, this is a non-trivial result. 

Because of the symmetry found for the MEP, a contour plot of the potential surface characterizes the 
path quite well. Constraining the H — H axis to be parallel to the surface and aligned with a surface groove, 
and fixing the Hi center of mass above a long-bridge site, the two remaining variables are Z, the H2 dimer's 
height above the surface, and R, the H-H separation. Figure 5 shows a contour plot of the H2/CU potential 
energy surface, V(Z, R), with the H2 geometry restricted in this way. The surface configuration is fixed at 
the minimum energy for a bare surface. The energies are in eV relative to the bare surface and two separated 
hydrogen atoms. 



6.4 Statistical sampling of the hyperplanes 

The hyperplane constraints derived from the MEP become quite simple due to the symmetry. Taking 
the x axis to be along the rows of atoms on the (110) surface, the y axis across the rows, and the z axis 
normal to the surface, the hyperplane constraint derived from the MEP becomes i\ + z 2 = constant in the 
reactant state and x\ — £2 = constant in the product state. Here x\ is the x coordinate of hydrogen l's 
centroid, etc. In the intermediate hyperplanes the constraint is 

a(xi - x 2 ) + b{zi + z 2 ) = 0. (40) 

where a and b are constants. This means that within the transition state hyperplane, for example, the 
centroids are allowed to move closer together only if they simultaneously descend toward the surface. They 
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are allowed to move apart only if they ascend higher above the surface. This constraint prevents the 
atoms from either chemisorbing on adjacent lattice sites, or desorbing. The plane constraint does not affect 
symmetric motion along the groove (a change in x\ + £2), asymmetric vertical motion (z\ — z-i) or the y 
coordinates perpendicular to the groove. The plane also does not prevent individual images from desorbing 
or adsorbing as much as the springs between them will allow. 

Figure 6 shows scatterplots of the H — H orientations in the classical statistical sampling at T — 100 K . 
Plotted are the angles theta and phi of the H-H orientation, relative to the surface. Theta is the polar angle, 
between and tt. Phi is the azimuthal angle, between — n and tt. Figure 6a is for a hyperplane in the 
reactant region. There is no restriction on the orientation of the molecule from the hyperplane constraint. 
There are fewer points at low and high values of theta because most of the surface area of a sphere is found 
near the equator, at theta « 7r/2. The scatterplot, when drawn as a rectangle, is therefore not uniform. In 
figure 6b is taken at the the transition state hyperplane where the rotation is highly constricted by the EAM 
interaction potential. This orientational confinement of the molecule in the transition state contributes to a 
significant entropy barrier for the transition. 

In the quantum statistical sampling, the two H atoms and eight Cu atoms nearest to the H atoms are 
included as FPI chains. We have used P = 50 images in each FPI chain in the present calculations. Figure 
7 shows three snapshots of the FPI chains from the quantum statistical sampling of the transition state 
hyperplane at T = 100 K. The system consists of two hydrogen atoms in the same groove of the Cu(110) 
surface, across the long bridge. The great flexibility available to the images in the H atom FPIs is evident. 
For illustration purposes only 20 beads where included in the FPI chains in this figure. 

The hyperplane constraint is a very loose constraint, allowing the system to sample a wide range of 
configurations. The system slides back and forth in the groove as illustrated in fig. 7. At low temperatures, 
where tunneling is important, the images drape down significantly on either side of the barrier, as illustrated 
by the separation of the images corresponding to the same FPI into two groups in fig 7. Few images are 
actually at the barrier top. At a higher temperature, the springs arc tighter and the distribution of images 
remains nearly spherical even at the barrier top. 



6. 5 The free energy barriers 

The accumulated free energy difference with respect to the reactant state as the hyperplane is moved 
along the reaction path in a classical statistical sampling is shown in fig. 8. The free energy rises faster than 
the potential energy along the reaction path. The free energy barrier increases with temperature, indicating 
a drop in entropy as the system approaches the transition state. The intercept of an Arrhenius line from 
several temperatures gives AS = —3.4 ks- 

This drop in entropy at the transition state is a result of the confinement of the molecule (see fig. 6) 
by the interaction with the surface as well as confinement of the lattice vibrations due to the presence of 
the molecule. The two effects were studied separately by first holding the lattice rigid and sampling only 
the degrees of freedom of the hydrogen molecule (fig. 9a); and then holding the hydrogen molecule fixed 
on the MEP while sampling statistically over the lattice degrees of freedom (fig. 9b). In the second case 
the rotational-vibrational entropy in the 5 in-plane degrees of freedom of the hydrogen molecule has been 
removed. The free energy becomes much closer to the MEP potential and has little temperature dependence, 
corresponding to AS = —1.2 fc^. The effect of the molecular confinement with the lattice rigid is much larger, 
corrsponding to AS = —5.3 ks- The two effects are not additive. The dynamics and relaxation (amounting 
to ca. 0.15 A displacement of the surface Cu atoms) reduce the confinement of the H 2 molecule in the 
transition state, thereby lowering the net entropy barrier for the dissociative sticking. In other words, the 
movable lattice restricts the hydrogen dimer by a smaller amount than does the rigid lattice. The lattice 
dynamics therefore assist the reaction by lowering the cntropic barrier. 

The results of the quantum statistical sampling are shown in Figure 10 for several different temperatures. 
The low temperature curves are flat-topped near the transition state. This is a tunnelling effect; it is a 
consequence of the chain's derealization. The images drape down over the barrier and manage to avoid the 
barrier top. Figure 11 compares the free energy curves for quantum H2, quantum D2, and the classical i?2, 
at 100 K. The deuterium curve shows clear tunneling effects, but smaller than hydrogen, as expected. The 
accumulated free energy in the quantum statistical sampling also drops below the classical free energy even 
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before the tunnelling region is reached. This is evidently a zero-point energy effect due to the softening of 
the H — H vibration. 

We have used a total of 45 hyperplanes to evaluate the free energy barrier. In the first hyperplane the 
center of the molecule is 7 A above the surface. The spacing between the hyperplanes was not even as can 
be seen from the data points in fig. 8. Initially the spacing is large since the force is small far from the Cu 
surface. The accuracy of this discretization of the hyperplane progression and the numerical integration of 
the force was tested by recalculating the free energy barrier using only every other hyperplane. At 100 K, 
the integrated free energy barrier differed by 0.0002 eV in the classical simulation and by 0.016 eV in the 
quantum simulation. The error due to finite hyperplane spacing is largest in the low-temperature quantum 
calculation where tunnelling is most significant, because the free energy curve is flat-topped and has a rapidly 
changing derivative at the onset of tunnelling. 

A larger source of error in the free energy calculation are statistical errors due to the limited simulation 
time in calculating the average force. This error is greatest in the high-temperature quantum simulations 
where the timestep needs to be very small. From the standard deviation in the force acting on the hyper- 
planes, we estimate the standard deviation in the accumulated value of the free energy barrier to be 0.03 eV 
in the quantum results at 600 K - which is the case with largest statistical fluctuations. 



6. 6 Tunneling 

In an effort to separate the tunneling effect from the effect of changes in zero-point energy, the RMS 
derealization perpendicular to the hyperplane was evaluated. The RMS derealization of a free atom in one 
dimension, as determined by recursive Gaussian integration 8 of the quantum partition function and taking 
the limit as P — > oo, is 

v/((x-x) 2 ) = h/y/l2mk B T, (41) 

i.e., inversely proportional to the square root of both temperature and mass. For a hydrogen atom at 100 
K, this is 0.2A. 

Figure 12a shows the average RMS derealization perpendicular to the hyperplane as a function of the 
reaction coordinate for H 2 at 100 K and at 300 K. At 100 K the derealization suddenly increases, in the 
space of one hyperplane, to a much higher value in the barrier region. This derealization is the source of 
the flat-topped free energy barriers at low temperatures and is a clear signature of tunneling. 

The temperature where tunneling becomes important can be estimated from the curvature of the po- 
tential barrier. Consider a chain with its centroid atop a one dimensional parabolic barrier. The images 
will spread away from the barrier top as far as the spring constant allows. The barrier top is thus not as 
unfavorable for a quantum atom as it is for a classical atom. The centroid density at the barrier top is 
correspondingly increased, and any reaction that goes over the barrier will proceed faster than expected. 
This corresponds to tunneling. 

The effective energy of a given configuration x(r), which determines the Boltzmann sampling of config- 
urations with V(x) a parabolic barrier potential V = —kx 2 /2 is 



Veff[x(r)]=J2 



>^spr\^i •^i—l) J- >v\3/i 



p 



(42) 



Both the spring energy and the barrier potential scale as distance squared. Thus, if the distances of all beads 
from the center are doubled (xi — » 2a;,), both the spring energy and the barrier energy are quadrupled. If 
the spring constant k spr is small compared to the barrier sharpness k, V e ff is negative. Then V e ff decreases 
without bound as the beads move away from the barrier, and the chain delocalizes to infinity in both 
directions. 

In practice, when the springs are weak, the chain runs away from the barrier top until the barrier starts 
curving back up. The critical temperature, in terms of the imaginary barrier frequency $7 = fc//i, turns out 
to be ksTc = Ml / (2ir) } 2 This is a reasonable estimate of the temperature at the onset of tunneling. Below 
T c , the chain derealization near the barrier is expected to be significantly larger than in free space. 

Fig. 12b shows the difference in RMS derealization between the reactant and transition state at a 
range of temperatures, for H2 and Z?2- Significant increase in the derealization sets in at about 400 K for 
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hydrogen and at about 300 K for deuterium, as the temperature dependent spring constants in the FPI 
become weaker than the curvature of the MEP potential barrier. Using the curvature of the MEP potential 
barrier to estimated the onset of increased derealization and the equation for T c above, gives T c — 400 K 
for H 2l and 290 K for D 2 , in good agreement with the observed behaviour. 

The onset temperature of tunneling derealization is well predicted from a parabolic barrier reference, 
but the derealization turns out to increases smoothly as the temperature is lowered further, while the 
parabolic barrier model predicts a very sharp cutoff temperature below which the derealization goes to 
infinity. The gradual increase in the derealization is due to the anharmonicity of the MEP potential. As the 
FPI images start probing regions where the barrier becomes significantly non-parabolic a balance is reached 
between the downward push from the potential gradient and the upward pull from the springs connecting to 
images on the other side of the barrier top. Thereby the tunneling derealization saturates. Fig. 13 shows 
the distribution of FPI images for hydrogen projected onto the MEP in typeical snapshots taken from the 
statistical sampling at 100 K and 300 K. Even at 300 K where the tunneling effect is small, the FPI images 
already probe regions where anharmonic effects are clearly significant. 

In H 2 at 300 K, the change in RMS derealization between initial and transition states (sec figure 
12b) is small. However, the actual derealization is about 0.12A The derealization goes to zero at high 
temperature, but very reluctantly, as an inverse square root. The quantum system therefore is not expected 
to become truly classical except at extremely high temperatures. 

However, the quantum system at moderate temperatures such as 300 K behaves quasi-classically: the 
derealization remains roughly constant, and the distribution of images remains roughly spherical. The 
system's statistics correspond roughly to those of a classical system at the same temperature, but whose 
potential energy surface has been shifted to account for the zero-point energy effects seen at the non-zero 
derealization. In other words, the quantum effects at moderate temperatures, though non-zero, may simply 
show up as an implicit shift in the potential energy, rather than as the fundamental differences that are 
observed at low temperatures. 



6. 7 Torque contribution to the reversible work 

The rotational contribution to the reversible work turns out to be small for this system. Figure 14 
shows the accumulated free energy difference due to rotation of the hypcrplanes as a function of the reaction 
coordinate. While the curvature of the path is quite large, the torque with respect to the MEP is small and 
the net effect on the free energy is small. This is unlike the Eckart barrier plus harmonic oscillator model 
system, where the rotational contribution was found to be just as large as the translational contribution 
when using the MEP as reaction path. 10 There, tunneling caused the average density to move off the MEP 
and thereby create torque on the hyperplane in a region of large path curvature. By redefining the path 
to coincide roughly with the average density, the torque on the hyperplane and thereby the rotational 
contribution to the free energy was minimized. The net free energy barrier, however, was the same. 

While the rotational contribution is small here, it is statistically significant, and of opposite sign in the 
classical and quantum case (fig. 14). Quantum derealization effects result in a shift of the average density 
of the system with respect to the MEP, which presumably changes the direction of the torque and thereby 
the sign of the rotational contribution to the free energy 



6.8 Activation energy for adsorption and desorption 

Fig. 15 shows an Arrhenius plot for the calculated free energy barriers to dissociative adsorption. The 
classical statistical sampling gives a straight line corresponding to a constant activation energy of 0.73 eV. 
Both H 2 and D 2 are observed to be classical at 600 K but below 400 K the quantum results deviate from 
the classical results. The quantum statistical sampling gives significantly lower free energy barriers below 
400 K, and a correspondingly lower activation energy of 0.38 eV for H 2 and 0.45 eV for D 2 . A near Arrhenius 
behaviour is observed in the quantum regime as well as in the classical regime; the calculated points fall 
close to a straight line fit. The lowering of the activation energy is mainly due to tunneling, as illustrated 
by the derealization perpendicular to the hyperplane shown in fig. 12b. The tunneling isotope effect is 
(0.73 -0.38)/(0.73- 0.45) = 1.25. 
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The experimental data of Campbell and coworkers is also shown in fig. 15. The measured activation 
energy agrees well with the calculated results, but the measured sticking probability is larger by a factor of 
40, indicating the entropy barrier is too large in the simulation and the model potential should be adjusted to 
have a wider reaction channel at the transition state. If we assume the in-plane vibrations are harmonic, this 
can be accomplished, for example, by reducing all five frequencies by a factor of 1/40 1 / 5 w 1/2. This takes 
into account only the five non-reactive hydrogen vibrations (such as symmetric translation in the groove 
of the (110) surface, symmetric and antisymmetric vibrations perpendicular to the groove, and antisym- 
metric vertical motion) and neglects modifications in the copper vibrations. Since the shape of our model 
potential surface is quite arbitrarily derived from the LEPS potential form with no direct experimental input 
determining the width of the reaction channel at the transition state, such a modification is not unreasonable. 

Fig. 16 shows an Arrhcnius plot for the reverse reaction, recombinative desorption. The classical 
statistical sampling gives Arrhcnius type behaviour with the calculated points falling near a straight line 
corresponding to a constant activation energy of 0.55 eV. The quantum statistical sampling gives signifi- 
cantly lower free energy barriers below 400 K, and a correspondingly lower effective activation energy. The 
calculated points do not fall on a straight line, as is typical for decay of metastable states, and the effec- 
tive activation energy is therefore somewhat temperature dependent. At 300 K the calculated activation 
energy is 0.49 eV for Hi- An experimental activation energy has been reported by Anger, Winkler, and 
Rcndulic 34 as 0.43 eV in this temperature range and in the zero coverage limit. The calculated and measured 
activation energies are therefore in quite good agreement. The activation energy at 100 K is calculated to 
be significantly lower, down to 0.36 eV. Again, tunnelling through the barrier has a strong effect at low 
temperature. 



6.8 Classical dynamical corrections 

The transmission coefficient, k, can be calculated by molecular dynamics. 8 ' 16 An ensemble of points 
is found at the dividing surface whose momenta are Boltzmann-distributed and initially forward-directed. 
The ensemble is propagated forward and backward in time until all its members leave the vicinity of the 
transition state. The fraction of these trajectories which actually started in R and finish in P is counted. 
We assume that once the system leaves the vicinity of the transition state, it will remain either a reactant 
or a product for a much longer time than it spent near the transition state. If, however, it spends a long 
time in the intermediate region, dynamical corrections become problematic. 

Classical dynamical correction calculations were carried out for H 2 and D 2 at 600 K, where classical 
statistical sampling gave the same TST rate estimate as quantum statistical sampling. An ensemble of 144 
points was taken at the transition state hyperplane; each of these was assigned a Boltzmann velocity and 
each trajectory was followed forward and backward in time until it left the transition state region. TST 
assumes that each crossing point will finish as products when propagated forward and as reactants when 
propagated backward. If, for a given crossing point, cither of these assumptions is untrue, it does not count 
toward the rate. The transmission coefficient was found to be k — 0.54 ± 0.08 for H 2 and k — 0.56 ± 0.08 
for D 2 . No kinetic isotope effect was observed for this system. 

A transmission coefficient of 0.5 has the same effect on the rate as does a free energy increase of 
—kBTlnn = 0.03 eV at 600 K. This is of the same magnitude as the error associated with calculation of 
the free energy. More careful estimations of k are therefore not called for unless the free energy difference is 
also calculated significantly more accurately. 

None of the trajectories lingered in the transition region for more than a picosecond. The quantum 
dynamics are expected to be similar to the classical dynamics at this high temperature, and the quantum 
transmission coefficient will, therefore, likely be within an order of magnitude of the classical transmission 
coefficient, corresponding to a few hundredths of an eV. The free energy difference would then provide 
a good estimate of the sticking coefficient, even without dynamical corrections. Furthermore, since the 
classical trajectories spent little time near the transition state, it may be feasible to do a quantum-dynamical 
calculation for the short time required to estimate corrections to the QTST estimate. However, we expect 
those corrections to be small, and conclude from the classical dynamical corrections that TST is an excellent 
approximation for treating the H 2 dissociative adsorption as well as associative desorption at a Cu(110) 
surface. Other dissociative adsorption processes are likely to be also well approximated by TST. 
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7. Conclusions 

We have presented a calculation of the dissociative sticking and thermal desorption of hydrogen molecules 
using a statistical approach, quantum transition state theory, where the two hydrogen atoms were included 
fully quantum mechanically and a thermal average was carried out over both the molecular and surface 
degrees of freedom. For the model potential surface used here, the system crosses over from classical to 
quantum behavior at about T — 400 K for Hi and T = 300 K for Z?2- In the quantum regime, the 
activation energy is significantly reduced, the free energy barriers are flat-topped, and the derealization is 
significantly greater at the barrier than in the reactant state, indicating that tunneling becomes the dominant 
transition mechanism. 

A practical method for evaluating the required free energy barriers in classical and quantum systems 
has also been given. This method is based on a reversible work evaluation and can be applied with relative 
ease to systems of high dimensionality. Another advantage of this method is the identification of an optimal, 
hyperplanar dividing surface for the TST rate estimate in the course of the reversible work calculation. 



17 



APPENDIX 

In this appendix we show that the contributions to the free energy difference between two hyperplancs 
due to the in plane forces vanishes in the quantum mechanical case, just as in the classical case. We start 
with the rotational component of cqn. (25) 



dArot (s) 



r^ h dr \ 

J ^VV[r(T)]-U' s Uj(r(r)-r s )\ da. 



(Al) 



With the transformation to local coordinates, z, given by eqn. (12), this expression may be rewritten in 
terms of the rotated coordinates z as 



f3h dT dV[z(T)} 

I, (3h dz 



dA rot (s) = I I ^' i :y n ■ UfU' s z(r) ) ds 



(A2) 



where the average is written in terms of z as 



(■■■>« 



/Dz(r)e- s / R ---d[(io)i] 



/£)z(r)e- s /^[(z )i] 
Consider the following average where i ^ 1 and j =/= i 

'■ f3n dr dV[z(r)} \ 



J-ij — 



(3h dzi 



- z j( T ) 



(A3) 



(A4) 



In terms of a path integral, this becomes 



J-ij — 



c J w* f Z'omwv.H 



(A5) 



where C is a normalization constant. 
Using the following relation 



dS 

d(z Q ) 



0h d _ dV(z(r)) 



we may write 



Ui = ~j\ IMr) 



d(z ) 



e- s / R ^(r)5[(I )i; 



Iij can be shown to be zero in the following way. First introduce 



(A6) 



G 



e- s / R ^(r)5[(z )i] 



(A7) 



so that In can be written as 



and 



— G f f) 

(3 J 'd{z )i 



G — > as (zo)i — ► ±oo. 



By transforming to Fourier variables 



(A8) 

(A9) 



1 



hrl 



z^-j Q dre-^z(r) 



(AW) 
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where 



eqn. (A6) becomes 



fit 



2-rrk 
JifJ 



—C 

ho = —J 



\{ Jd{i k )i fd(z 



(fe,0^(o,» 
where J is a constant Jacobian for the transformation 

f Dz(t) = J H J d(z fc ),. 

We next do the (zq)j integration to give 



d(z ) 



G 



(k,i)- 



—C 
hi = —J 



n /<w 



(k,l)jt(0,i) 



(G(+oo) - G(-oo)) ds 



which verifies 



^=0. 



Using this result and the fact that Ug U' s is antisymmetric, gives 



dA rot (s) 
which can be rewritten as 

dA ro t(s) = - 



n^M *° 



? h dT_ 

o /3?i 



g^[r(r)] 
Or 



•U, 



[uflrM-r,^ 



and with [U s ] i:L = [n s ] i7 we recover our final result 

iA t s / f m dr f dV[r(r)} 



n s « • (r(r) - T s )) ) ds 



which is equivalent to eqn (27) derived for the classical case. 



(All) 
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potential and the coupling to the thermal bath between the present calculations and ref. 9 leading to 
small differences in the calculated results. 

Table 1: Properties of the empirical interaction potential 



Characteristic 

Separation (A) 
Binding energy (eV) 
Vibration (eV/A 2 ) 



Binding energy (eV) 7,34 
Vibration (eV/A 2 ) 38 

Saddle point energy (eV) 
H-H separation (A) 
H-Cu height (A) 



Experimental 

0.737 

-4.747 

24.50 



Ho 



H-cu(noy. 



-2.3 

5.5 



Long Bridge Saddle Point: 



Fitted 

0.740 

-4.7472 

25.43 



-2.29 
5.0 



+0.688 
1.3385 
0.8802 
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Cu 


H 


2862 


79.50 


3.512 


2.480 


-109.1 


-107.6 


1.756 


2.999 


0.273 


2.144 


3.691 


3.777 


7.381 





6 





512 





-112.9 


-81.75 


8510 


838.7 


-2.617 E+05 


-3953 


4.780 E+06 


8768 


-5.234 E+07 


-6599 


3.391 E+08 




-1.201 E+09 




1.796 E+09 





Table 2: Parameters of the interaction potential (units of eV and A) 

Parameter Cu H Cu-H 

D A 2862 79.50 86.15 

a A 3.512 2.480 4.211 

D B -109.1 -107.6 1.536 El 04 

a B 1.756 2.999 6.076 

S 

Pa 

8b 

') 

7 

h 

h 

h 

h 
h 
h 

h 
h 

Figure captions 



Fig. 1. (a) Reaction path T parametrized with reaction coordinate s. For each point along the path, a 
hyperplane intersecting the path at that point and having a normal n s is defined, Z = n s • (r — r s ) = 0. 
The force acting on the system is averaged over a statistical sample of configurations where the system is 
constrained to lie in the hyperplane. 

(b) The hyperplane is gradually moved from the reactant region towards products by varying the reaction 
coordinate, s. The activation free energy for the transition is obtained by calculating the reversible work 
involved in moving and rotating the hyperplane to an intermediate position chosen to be the dividing surface 
between reactant and product regions. The maximum free energy barrier is obtained where the force and 
torque on the hyperplane vanish. This is equivalent to finding the optimal hypcrplanar dividing surface in 
classical variational transition state theory. 



Fig. 2. The system confined to a hyperplane experiences a gradual change in the potential energy as the 
hyperplane is moved along the reaction path. (a) A schematic contour plot of a potential energy surface 
with the reaction path shown as a thick solid line. Two different locations of the hyperplane are indicated 
along the progression from reactants to products. They intersect the reaction path at T s and r s+( 2s and 
their orientation is given by the normal vectors n s and n s+ ds as shown in the inset, (b) Potential energy 
curves corresponding to the two hyperplanes. This illustrates the ramping of the potential experienced by 
the system confined to stay within the hyperplane, as the hyperplane is moved from one location to another. 



Fig. 3. In the quantum mechanical evaluation of the free energy barrier the hyperplane constraint is 
applied to the centroids, ?o, of the Fcynman path integral representations of the quantum particles in the 
system, instead of the classical atom coordinates (compare with fig. la). The force acting on the FPI images 
is evaluated and averaged over a statistical sample of configurations, with the centroids of the FPI chains 
constrained to lie in the hyperplane. 



Fig. 4. The minimum energy path between the molecular hydrogen outside the surface and the dissociated 
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hydrogen chemisorbed on adjacent sites on the Cu(llO) surface. The Nudged Elastic Band method was used 
to find the optimal path. The initial state, shown in the upper left panel, was arbitrarily chosen to have 
the molecule oriented with its axis normal to the surface plane and located directly above a surface atom. 
The final state, shown in the upper right panel, has the H atoms sitting in fourfold sites separated by a 
long bridge. Along the minimum energy path, shown in the lower panel, the molecule shifts over to the the 
long bridge site and turns to make the axis parallel to the surface. The path from then on is symmetric and 
involves motion perpendicular to the surface and perpendicular to the long bridge. The molecule breaks up 
very close to the surface. 



Fig. 5. Potential energy contours for the dissociation of H 2 on Cu(llO) over a long bridge site. Because of 
the high symmetry of the MEP shown in fig. 4, a two dimensional cut through the six dimensional surface 
effectively describes the energetics. The molecule is coming down parallel to the surface and is symmetrically 
arranged over a long bridge site; z is the height above the surface and R the intramolecular distance. 



Fig. 6. Scatter plot of the orientation of the H^ molecule in the statistical sampling of the classical system 
at T = 100 K. Theta is the polar angle, between and it. Phi is the azimuthal angle, between —it and n. 
(a) Points taken from the sampling of a hyperplane near the reactant region. The hyperplane constraint 
does not restrict the rotation of the molecule, (b) Points taken from the sampling of the transition state 
hyperplane, illustrating the rotational confinement of the molecule due to the Hi — Cu interaction potential. 



Fig. 7. Snapshots from the statistical averaging of the transition state hyperplane at T = 100 K. 
Only one of the six centroid degrees of freedom of the H atoms is constrained, so the system still has five 
degrees of freedom to explore the potential energy surface. In particular, the hyperplane constraint does 
not prevent simultaneous movement of the H atoms perpendicular to the bridge, as is evident from the 
snapshot shown in the top panel. However, the centroids are not allowed to move closer together unless they 
simultaneously descend toward the surface. They are allowed to move apart only if they ascend higher above 
the surface. This hyperplane constraint prevents the atoms from either chemisorbing on adjacent lattice 
sites, or desorbing. The images of each of the FPI chains tend to separate into two clusters illustrating the 
sliding down on opposite sides of the potential barrier. This is a signature of tunneling. 



Fig. 8. The accumulated free energy difference with respect to the reactant state as the hyperplane 
is moved along the reaction path in a classical statistical sampling. For reference the solid line shows the 
potential energy along the minimum energy path. The free energy barrier increases with temperature, 
indicating a decrease in entropy in the transition state. 



Fig. 9. Free energy curves obtained in classical statistical sampling after freezing out selected degrees 
of freedom of the system. (a) The lattice is made rigid to illustrate the effect of confinement of the Hi 
molecule as the transition state is approached. The entropy barrier becomes larger when the lattice is rigid 
(compare with fig. 8), illustrating that relaxation in the lattice relieves the confinement of the molecule in the 
transition state, (b) The dimer is held rigid on the reaction path while the lattice degrees of freedom are 
sampled. The molecule confines the lattice vibrations significantly only at high temperatures, T > 300 K . 



Fig. 10. Quantum statistical sampling of the hyperplane along the reaction path (quantum mechanical 
extension of the results in fig. 8) showing a significantly lower free energy barrier at low temperature as 
compared with the classical statistical sampling. Already at T — 300 K quantum effects have lowered the 
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free energy barrier by ca. 0.1 eV. At 100 K the barrier is flat-topped due to large tunneling effects. The free 
energy curve drops below the potential energy curve early along the reaction path, well before tunneling is 
effecctive, due to lowering of the zero point energy of the molecule. 



Fig. 11. The accumulated free energy at T = 100 K for a quantum statistical sampling of D 2 dissociative 
adsorption The classical and quantum sampling of H 2 dissociative adsorption from figs. 8 and 10 is shown 
for comparison. The classical sampling gives the same results for D 2 and H 2 . This shows that significant 
lowering of the free energy barrier is found for D 2 , but smaller than for H 2 . 



Fig. 12. (a) The RMS derealization of H 2 perpendicular to the hyperplane as a function of the reaction 
coordinate. In a free H atom the derealization is 0.2 A at 100 K . As the repulsive interaction between the 
molecule and the surface sets in the derealization decreases. Near the transition state, a dramatic increase in 
the derealization is observed at 100 K due to tunneling. In this segment of the reaction path, the free energy 
is nearly constant (see fig. 11). At 300 K, the derealization is much weaker but still noticeable, (b) The 
difference in the RMS derealization perpendicular to the transition state hyperplane and the derealization 
of free H 2 molecule and D 2 molecule as a function of inverse temperature. Straight lines fitted to the larger 
values are shown only to guide the eye. Significant increase in the derealization sets in at about 400 K 
for hydrogen and at about 300 K for deuterium, as the temperature dependent spring constants in the FPI 
become weaker than the curvature of the MEP potential barrier. 



Fig. 13. Distribution of FPI images for hydrogen projected onto the MEP in snapshots taken from the 
statistical sampling at 100 K and 300 K. The dashed line shows a parabolic barrier with curvature matching 
that of the MEP potential curve. Even at 300 K the FPI images probe regions where anharmonic effects 
are significant. The vertical position of the images is arbitrary. At 300 K the energy of the system is high, 
several eV over the barrier, and the distribution of images does not reflect the shape of the barrier. At 100 K 
the energy of the system is closer to the barrier energy and the effect of the barrier on the distribution is 
clear. 



Fig. 14. Rotational contribution to the free energy as a function of reaction coordinate at T = 100 K. 
The contribution is negligible as compared with the translational contribution in this system when the the 
MEP is chosen to be the reaction path, because the torque acting on the hyperplanes turns out to be very 
small. A small but statistically significant contribution is obtained where the path curvature is largest, near 
the transition state. Quantum derealization effects result in a shift of the average density of the system 
with respect to the MEP which presumably changes the direction of the torque and thereby the sign of the 
rotational contribution to the free energy. 



Fig. 15. An Arrhenius plot for H 2 and D 2 dissociative adsorption. The calculated points are for 100 K, 
150 K, 200 K, 300 K, 400 K, and 600 K. The classical statistical sampling gives a straight line corresponding 
to a constant activation energy of 0.73 eV. The quantum statistical sampling gives significantly lower free 
energy barriers below 400 K and a correspondingly lower activation energy of 0.38 eV for H 2 (slope of 
large dashed line) and 0.45 eV for D 2 (slope of small dashed line). The lowering of the activation energy is 
mainly due to tunneling. The experimental data of Campbell and coworkers 7 is shown with triangles. The 
measured activation energy agrees well with the calculated results, but the measured sticking probability is 
larger by a factor of 40, indicating the entropy barrier is too large in the simulation and the model potential 
should be adjusted to have a wider reaction channel at the transition state. A solid line is drawn through 
the calculated points only to guide the eye. 
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Fig. 16. An Arrhenius plot for H2 and D2 associative desorption. The classical statistical sampling gives 
a straight line corresponding to a constant activation energy of 0.55 eV. The quantum statistical sampling 
gives significantly lower free energy barriers below 400 K , and a correspondingly lower activation energy. At 
300 K, the calculated activation energy is 0.49 eV for H2 as compared to the experimental value 0.43 eV 
obtained by Anger et al. (the solid line has a slope consistent with the experimental measurements but its 
vertical position was simply chosen in such a way as to intersect the origin, since an absolute value of the 
desorption rate was not reported). A straight dashed line is drawn through the calculated points only to 
guide the eye. 
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